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Understanding the detailed production and hadronization mechanisms for heavy quarkonia and 
their modification in a nuclear environment presents one of the major challenges in QCD. Calcula- 
tions including nuclear-modified parton distribution functions (nPDFs) and fitting of break-up cross 
sections (atr) as parameters have been successful at describing many features of J/ip modifications 
in proton(deuteron)-nucleus collisions. In this paper, we extend these calculations to explore differ- 
ent geometric dependencies of the modifications and confront them with new experimental results 
from the PHENIX experiment. We find that no combination of nPDFs and atr, regardless of the 
nPDF parameter set and the assumed geometric dependence, can simultaneously describe the entire 
rapidity and centrality dependence of J/ip modifications in d+Au collisions at <Js NN = 200 GeV. 
We also compare the data with coherence calculations and find them unable to describe the full 
rapidity and centrality dependence as well. We discuss how these calculations might be extended 
and further tested, in addition to discussing other physics mechanisms including initial-state parton 
energy loss. 



I. INTRODUCTION 

In proton(deuteron)-nucleus collisions at the Relativis- 
tic Heavy Ion Collider (RHIC), the nucleus is extremely 
Lorentz-contracted and thus the entire interaction and 
traversal of the nuclear target takes place on a time scale 
of order 0.1 fm/c. Thus, one expects coherence effects to 
play a significant role in the physics of particle produc- 
tion and hadronization. By studying heavy quarkonia 
states, one can postulate that the initial hard produc- 
tion of a cc pair can be factorized from the later traver- 
sal of that pair through the remainder of the nucleus. 
Many calculations have utilized this factorized approach 
in trying to understand the nuclear modification of J/ip 
yields in proton(deuteron)-nucleus reactions (for exam- 
ple [IH1]). In this factorized framework, modification of 
initial cc pair production is accounted for via nuclear- 
modified parton distribution functions (nPDFs). After 
being produced, the disassociation of charm pairs and 
thus the additional reduction in the number of final-state 
J/ip mesons is accounted for with a simple breakup cross 
section (<7(, r ). It is interesting to test this picture to see 
if the experimental data requires additional physics, in- 
cluding coherence effects and initial-state parton energy 
loss. This work presents an extension of this simple cal- 
culational framework and investigation of some of the key 
underlying assumptions. 
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II. CALCULATION DETAILS 

In this section we describe the inputs required for the 
calculation of the nuclear modification factors (R p a or 
RdA) for various nuclear targets and centrality selections. 
First, the density of partons in the nucleus is modified 
relative to the parton distribution function (PDF) for 
nucleons. This nuclear-modified PDF reflects the modi- 
fied parton density encountered by the projectile and re- 
sults in a modified number of hard scatterings that create 
cc pairs from g + g, q + g, and q + q interactions. The 
state of the art calculation is the EPS09 nPDF parameter 
set with uncertainties represented by 31 different Hessian 
basis parameterizations as detailed in [Sj. Because J/ip 
production at high energies is dominated by interactions 
between gluons, we will consider g+g interactions only in 
the calculations in this paper. Figure [T] shows the EPS09 
gluon modification Rq at a Q 2 = 9 GeV 2 , the appro- 
priate scale for production of the J/ip. It can be seen 
that the nPDFs are not well constrained by experimen- 
tal data, particularly the low-x gluon distributions which 
dominate the J/ip production probability at forward ra- 
pidity at RHIC energies. 

The second main effect is that after creating the cc pair 
in the initial state (often referred to as the J/ip precursor, 
since the hadronic state is expected to take of order 0.3 
fm/c to form), the pair may break up or be de-correlated 
while traversing the remaining portion of the nucleus. 
This second effect is often included by assuming a fixed 
cross section a^r for the breakup of the pair. We note 
that this effect is also often termed absorption, though 
this nomenclature can be misleading since the charm pair 
still exists, but is no longer able to form a final-state J/ip 
meson. Currently there is no fundamental description of 
the hadronization process for the J/ip meson that agrees 
with all the available experimental data [7j . The lack 
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of such a theory for the dynamics of hadronization means 
one has no ab initio calculations of this precursor-nucleon 
cross section and its dependence on the relative velocity 
between the pair and the target nucleons. In most works, 
the value of Obr is assumed to be independent of the J/ip 
rapidity for a given y?s NN , and is determined from fits to 
the experimental data |5]- 
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FIG. 1. The gluon nuclear modification Rg for the Au nu- 
cleus at the scale Q 2 — 9 GeV 2 is shown for the EPS09 cen- 
tral value (labeled set 1) and for all 30 error sets. The yellow 
shaded area is the overall uncertainty band calculated from 
the error sets, representing a 90% CL uncertainty. The pa- 
rameter sets from top to bottom at the lowest x values are 
16, 10, 7, 8, 13, 31, 26, 18, 15, 4, 3, 21, 23, 25, 28, 1, 29, 2, 
24, 20, 22, 14, 5, 19, 27, 30, 6, 12, 9, 11, 17. 

We employ a Monte Carlo Glauber model [9 J where the 
nucleons are randomly given spatial distributions within 
the deuteron based on the Hulthen wave function and for 
the gold nucleus based on a Woods-Saxon distribution 
with parameters R=6.38 fm and a=0.54 fm [2]. Individ- 
ual d+An collisions at y/s NN = 200 GeV are simulated 
by randomly selecting an event impact parameter (6) and 
determining if any pair of nucleons collide using an in- 
elastic cross section of a = 42 mb. One example event is 
shown in Figure [2] where the open circles are the posi- 
tions of the gold nucleus nucleons in the transverse plane, 
the red filled circles are the positions of the two nucle- 
ons from the deuteron, and the green filled circles are the 
gold nucleus nucleons which suffered a binary collision. 
Each binary collision between a red circle (deuteron nu- 
cleon) and a green circle (gold nucleon) has a probability 
to produce a cc pair. This probability is modified from 
proton-proton collisions according to the afore-mentioned 
nPDFs. 



The EPS09 nPDF parameterization, as well as other 
nPDF parameterizations, are predominantly determined 
from deep inelastic scattering experiments and minimum 
bias p + A reactions producing Drell-Yan pairs [5]. In 
such experiments there is no measure of the impact pa- 
rameter or transverse distance within the nucleus for the 
interaction and therefore the geometric dependence of 
the nPDF modification is not constrained. One expects 
a significant geometric dependence in the nPDF modi- 
fication, with a stronger modification near the center of 
the nucleus where the density of nucleons is the largest. 
In a simple picture, the partons inside the nucleons at 
low- a; have wave functions that are longer in the longi- 
tudinal direction than the Lorentz contracted nucleus. 
Thus the nPDF modification depends on the density of 
overlapping nucleons as shown in the transverse plane in 
Figure [2j However, there is no such Lorentz contraction 
in the transverse direction, and the parton wavefunction 
extent in this plane is of order 1 fm. Therefore, the largest 
nuclear effect would be expected where the density, and 
thus the longitudinal overlap is largest, near the center of 
the nucleus and should decrease as one moves out toward 
the periphery. 
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FIG. 2. Monte Carlo Glauber event display in the transverse 
(x-y) plane. Open black circles are the positions of the gold 
nucleus nucleons. Red filled circles are the positions of the 
two nucleons from the deuteron. Green filled circles are the 
positions of the gold nucleus nucleons which suffer at least 
one binary collision. The dashed black circle represents gold 
nucleus extent from the Woods-Saxon parameterization. 

In [10] , the nPDF modification is postulated to be lin- 
early proportional to the density-weighted longitudinal 
thickness of the nucleus at the transverse position of the 
binary collision, as written below: 

M(r T ) = 1.0 - aA(r T ) (1) 

where A(r r ) = — J dzp(z,r T ) is the density-weighted 
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longitudinal thickness and po is the density at the center 
of the nucleus. In Figure [2j each green circle is a trans- 
verse distance r T from the center of the gold nucleus, and 
one can then determine the average local thickness A(r x ) 
based on the Woods-Saxon parameterization. 

As shown in Figure [2j there are significant fluctua- 
tions in the thickness A(r T ) due to the randomly se- 
lected spatial locations of the nucleons in the gold nu- 
cleus at the time of the collision. In fact, the inclusion 
of such fluctuations has proven to be crucial in model- 
ing the initial conditions in heavy ion reactions (see for 
example [TTMl5| y In order to incorporate these fluctua- 
tions, we calculate the number of target nucleons around 
the struck nucleon (green circle) within a transverse ra- 
dius Rtube — 2 X Rnucleoni 

where R nuc ieon is taken to 
be the charge radius of the proton, 0.87 fm [16 . The 
average number of nucleons within this cylinder N tu be 
is it Rtube" 2 Po times A(r T ), and this proportionality con- 
stant can be absorbed into the parameter a in Eqn. [T] 
The exact choice of Rtube is somewhat arbitrary, but rea- 
sonable changes in the value do not significantly change 
the results shown in this paper. 
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FIG. 3. Shown is the Monte Carlo Glauber result for the 
Ntube distribution as a function of r T . The mean values of 
Ntube as a function of r T are shown as white points. The 
analytic calculation of the average value of A(r T ) from the 
Woods-Saxon parameterization, rescaled to the N tu be value, 
is shown as the solid red curve. Smearing the analytic calcu- 
lation around r T by the tube radius Rtube is shown as the red 
dashed curve. 

In Figure [3] the resulting distribution of N tu be values 
as a function of r T is shown. For collisions in the middle 
of the nucleus (r T 0) (N tube ) 20 and the RMS 
5. Using the analytic calculation of the average A(r T ) 
from the Woods-Saxon parameterization times n Rtube' 'Po 
yields the solid red curve. The difference between the two 
is because the N tu t, e {r T ) calculation includes the density 
averaged over the tube radius Rtube- If we smear the 



analytic calculation around r T by the tube radius Rtube, 
we obtain the red dashed curve, which now shows much 
better agreement. 

In our calculations, we take these fluctuations into ac- 
count by utilizing N tu be instead of the average A(r T ), as 
shown in the equation below: 

M{r T ) = 1.0 -aN tube (r T ) (2) 

There is a direct relationship between the parameter a 
and the (M), the modification averaged over all nuclear 
geometries, which is equivalent to i?dAu(0-100%). This 
relationship is determined by averaging M over the r T 
distribution for unbiased collisions as determined from 
the Monte Carlo Glauber (shown in Figure 3 of the 
PHENIX publication [17] ). The results are shown as the 
red curve in Figure[4] We also consider two other geomet- 
ric dependencies for the nPDF, referred to as exponential 
and quadratic. 

M(r T ) = Exp{-aN t ube{r T )) (3) 
M(r T ) = 1.0-a[iV ftlbe (r T )] 2 (4) 

The same procedure is applied, and the dependence of 
the parameter a is also shown in Figure [4] as the black 
(blue) curve for the exponential (quadratic) case. 
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FIG. 4. Shown is the parameter a as a function of the 
-RdAu (0-100%) modification that is obtained. The red, black, 
and blue curves correspond to the linear, exponential, and 
quadratic geometric dependence cases respectively. 

We examine the modification M(r T ) for each of these 
cases. Shown in Figure [5] (left panel) is M(r T ) for the ex- 
ponential case. The four solid black curves correspond to 
geometry averaged modifications of (M) = 0.8, 0.6, 0.4, 
0.2 for the top to bottom curves. The dot-dashed black 
line is the shape of the mean of the N tu b e distribution 
versus r T , normalized to one at r T = for reference. The 
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FIG. 5. Shown are the modification dependencies M{r T ) for the three geometry cases (exponential, linear, and quadratic) 
assuming four different average RdAu (0-100%) values - 0.8, 0.6, 0.4, 0.2 (corresponding to the top to bottom curves in each 
case) . The dot-dashed black line is the shape of the Ntube distribution versus r T , normalized to one at r T = for reference. 



middle (right) panel shows the same results for the lin- 
ear (quadratic) case. As a consequence of the linear and 
quadratic functional forms, M(r T ) has negative values 
for small r T when the geometry-averaged modification is 
less than 0.4 (0.6) for the linear (quadratic) cases. This 
unphysical result can be removed by forcing the modifi- 
cation to be positive-definite, but then one has to recal- 
culate the corresponding a parameters. In the analysis 
presented in this paper, the modification values from the 
nPDF do not typically reach these low values and thus we 
have not had to recalculate the results. However careful 
attention to this problem will be crucial for cases with 
larger modifications (for example more forward rapidity 
J/ip measurements or very low-x measurements at the 
Large Hadron Collider in p + A and A + A). 

There are two final ingredients needed to map the 
nPDF modifications onto the final state J/ip suppres- 
sion. These are the distribution of Bjorken X2 and Q 2 
for the parton-parton processes that contribute to the 
J/tp production and the mixture of g + g, g + q, and 
q + q processes. The simplest relationship between the 
J ftp and partonic kinematics arises under the assump- 
tion that the production is a 2 — > 1 process, for example 
g + g — > J/ip. Invoking conservation of energy and mo- 
mentum the production of a J/ip with p T = GeV/c 
results in the following relationship between J/ip rapid- 
ity (y) and parton momentum fraction (a^)- 
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This 2 — > 1 process is actually forbidden by angular mo- 
mentum conservation, but may approximate the correct 
kinematics at low p T , or in a color evaporation picture 
where soft gluon emission does not significantly modify 
the exact correlation of X2 and y. It has been pointed out 
that with a more detailed understanding of the subpro- 
cesses that contribute to J/ip production, one can utilize 
a more exact map of Xi and Q 2 to the final J/ip as a 



function of rapidity and p T |18H20j . The authors utilize 
the following relation between x\ and X2 that requires a 
full modeling of the cross section dependencies: 



x 2 = -, (6) 

[V^ x ^ - Vpt 2 + M 2 e yj 

This relation is exact for a 2 —> 2 process where one 
outgoing particle is an on-shell J/tp and the other particle 
is massless. 

Shown in Figure [6] (upper panel) is the correlation of 
Bjorken X2 with the J/ip rapidity in the 2 — > 1 case. 
This is compared with a scatter-plot showing calculation 
results from PYTHIA 6.416 [21] with the NRQCD setting 
for J/ip production. As expected, the x 2 values for a 
given J/ip rapidity are shifted to larger values. Since 
the J/ip (p T ) w 2.2 GeV/c, there must be a balancing 
particle(s), which requires larger available energy. Also, 
the emission of a balancing gluon, for example, will smear 
the rapidity of the J/ip relative to the 2 — > 1 calculation. 
Also shown are the (#2) values as a function of rapidity 
for the PYTHIA g + g color singlet channel (black) , g + 
g color octet channel (blue), and the q + g color octet 
channel (red). The mean X2 value can be misleading 
since it may have a large influence from a small fraction 
of high- a; events. Thus, in the lower panel we show the 
log(x2) distribution for the J/ip rapidity 2.0 < y < 2.4 
for the three different contributions. The majority of 
processes for this rapidity involve X2 « 0.002, but with a 
more significant high-x tail in the octet cases. Note that 
the underlying PYTHIA production does not obey the 
2—^2 kinematics of Eqn. |6j since there is initial-state fc-r 
and many of the octet production channels involve more 
than two final-state particles. 

We now incorporate all of the following items: (1) 
Monte Carlo Glauber, (2) deuteron and gold nuclear ge- 
ometries, (3) EPS09 nPDF parameter set, (4) geometric 
dependence assumption for nPDF (i.e. linear, quadratic, 
exponential), (5) kinematic mapping (g, q, q, X2, Q 2 — > 
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FIG. 6. (Upper panel) The kinematic correlations between 
the J/ip rapidity and Bjorken x% for the 2 — > 1 process com- 
pared with a scatter-plot showing calculations from PYTHIA 
6.416 with the NRQCD setting for J/ip production. The (x 2 ) 
values are shown for three different production mechanisms 
within PYTHIA 6.416. (Lower panel) Shown are the x 2 dis- 
tributions for J/ip produced with rapidity 2.0 < y < 2.2 for 
the three different PYTHIA production mechanisms (black 
for g + g color singlet, blue for g + g color octet, and red for 
q + g color octet). 



J/ip y, p T ). We then add the second factorized part of 
the calculation for the a^. by checking all nucleons on 
the back side of the nucleus (i.e. z nuc i eon > zunary) for 
whether the cc pair breaks up. Note that although the 
calculations are factorized, the results are auto-correlated 
by the geometry. For example, a binary collision occur- 
ring near r T = has a larger nPDF modification and also 
a larger probability for breakup. These auto-correlations 
are important to account for, and have been previously 
explored in terms of kicks broadening the p T distri- 
bution mug. 



Two additional benefits of this Monte Carlo Glauber 
approach with full fluctuations are that we can model 
the exact PHENIX experimental d+Au centrality selec- 
tion event-by-event and that we never project onto an 
averaged quantity (eg. the average impact parameter for 
each centrality class) and then calculate the modification 
for that average quantity. 




JAj» Rapidity 

FIG. 7. Shown is the J/ip nuclear modification factor RdAu 
for 0-100% interaction centrality as a function of rapidity. 
The calculations include the EPS09 nPDFs with the linear 
geometric dependence and 2 — ► 1 kinematics. The yellow 
band shows the limits from all 31 EPS09 nPDF variations, 
but should not be interpreted as a one-standard deviation 
uncertainty band. The PHENIX experimental data are shown 
as black points. The lines are the point-to-point uncorrelated 
uncertainties and the boxes are the point-to-point correlated 
systematic uncertainties. Not shown is the additional ±7.8% 
global scale uncertainty. 



III. CALCULATION RESULTS 

Putting all these pieces together, we show an example 
calculation in Figure [7] of the J / ip nuclear modification 
RdAu a s a function of rapidity for 0-100% of the d+Au 
inelastic cross section. In this example, we utilize the 
2 — > 1 exact process mapping and the linear geometric 
dependence of the nPDFs. We show the default EPS09 
result with a fjf, r = 4 mb (red curve), all 30 other varia- 
tions for EPS09 with a e>v = 4 mb (black curves), and a 
calculation assuming no nPDF modification with = 4 
mb (blue curve). The PHENIX experiment has recently 
reported high statistics J/ip rf+Au at y/s NN — 200 GeV 
nuclear modification factors as a function of rapidity |17j 
which are also shown in Figure [7J Within systematic 
uncertainties from the experimental data and theoretical 
calculation, reasonable agreement is obtained. 
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Centrality 60-88% 
I Global Scale Uncertainty ±10.0% 




_ Centrality 0-20% 
Z Global Scale Uncertainty ±8.5% 




Centrality 0-20%/60-88% 
Global Scale Uncertainty ±8.2% 




FIG. 8. Comparison of PHENIX data with calculations (red 
curves) using the 31 EPS09 parameter sets with linear nuclear 
thickness dependence. Each EPS09 parameter set is shown for 
its own best fit <r or value. The blue curve shows the best fit 
for the best set, obtained with EPS09 set 17 and a^ r = 3.2 
mb. 



FIG. 9. Comparison of PHENIX data with calculations (red 
curves) using the 31 EPS09 parameter sets with quadratic 
nuclear thickness dependence. Each EPS09 parameter set is 
shown for its own best fit Obr value. The blue curve shows 
the best fit for the best set, obtained with EPS09 set 30 and 
Obr = 3.4 mb. 



We emphasize that this calculation utilized the 2 — > 1 
kinematics. We have performed the same calculation us- 
ing the various PYTHIA kinematics and find only a very 
modest decrease (smoothing out) of the rapidity depen- 
dence. The rapidity dependence in the calculation comes 
entirely from the nPDF dependence on Xi and Q 2 . Thus, 
utilizing the PYTHIA kinematics leads to a slight blur- 
ring of this relation and a general shift to larger X2 values, 
as expected from Figure [6] This flattening of the RdAu 
versus rapidity goes in the opposite direction of the ex- 
perimental data; however, the uncertainties in the nPDFs 
and other physics do not allow for any conclusion about 
the underlying production process. 

Shown in Figure [8] are the PHENIX experimental re- 
sults for RdAu for peripheral events 60-88% (top panel), 
^dAu for central events 0-20% (middle panel), and the 
ratio between them R^p 0-20% / 60-88% (lower panel). 
Note that the significant systematic uncertainties that 
can modify the rapidity dependence of the modification, 
referred to by PHENIX as type-B systematics, largely 
cancel in the i?cp ratio. 

We utilize the r T distributions for each centrality class 
shown in Figure 3 of the PHENIX publication [T7] to 
compute the expected modification in each centrality. 
There are many different statistical fits one can perform 



between the experimental data and our theoretical cal- 
culations. In this case, we perform a modified-^ 2 (x 2 ) fit 
to just the i?cp data (which provides by far the best con- 
straint on the rapidity dependence). The \ 2 fit method 
that accounts for both statistical and systematic uncer- 
tainties is detailed in [2"4"] . 

We have chosen again to utilize the 2 — > 1 kinematics 
and consider the linear, quadratic, and exponential geo- 
metric dependencies for the nPDFs. Shown in Figure [8] 
are the results for the linear geometric dependence. Each 
red curve represents one of the 31 EPS09 nPDF parame- 
ter sets and the best fit a br value for that parameter set 
(ie. minimum x 2 ) to the Rcp data. The calculation cor- 
responding to that red curve is then also shown for RdAu 
peripheral and central in the upper and middle panels, 
respectively. The blue curve represents the best fit over- 
all for all combinations for EPS09 parameter sets and a br 
values - corresponding to EPS09 nPDF parameter set 17 
and a a br = 3.2 mb. However, even this best fit has a x 2 
= 41.5 which corresponds to an extremely poor fit (i.e. 
probability less than 10 -4 ). The result with the geomet- 
ric nPDF exponential case are similar with the best fit 
from EPS09 nPDF parameter set 17 and a br = 4.2 mb 
and a poor % 2 = 50.6. 

Shown in Figure [9] are the results for the quadratic ge- 
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FIG. 10. The points are the PHENIX J/ip R C p versus RdAu- 
The ellipses are the one-standard deviation contours for the 
systematic uncertainties. The open circles, closed squares, 
and closed circles are from backward, mid, and forward ra- 
pidity respectively. The black curves are analytic calculations 
assuming a purely exponential, linear, and quadratic geomet- 
ric dependence for the nuclear modification 17 . The blue 
lines are the full calculation, using an exponential thickness 
dependence for the shadowing, for all EPS09 nPDF parame- 
ter sets with cjbr = 0. The red lines are the full calculation 
for all EPS09 nPDF parameter sets with all at r values. 



Linear nPDF Case: All EPS09 sets + Alia 
Linear nPDF Case: All EPS09 sets +o u =0.0 mb 
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FIG. 11. The points are the PHENIX J /if) R C p versus RdAu- 
The ellipses are the one-standard deviation contours for the 
systematic uncertainties. The open circles, closed squares, 
and closed circles are from backward, mid, and forward ra- 
pidity respectively. The black curves are analytic calculations 
assuming a purely exponential, linear, and quadratic geomet- 
ric dependence for the nuclear modification 17 . The blue 
lines are the full calculation, using a linear thickness depen- 
dence for the shadowing, for all EPS09 nPDF parameter sets 
with Gbr = 0. The red lines are the full calculation for all 
EPS09 nPDF parameter sets with all <7{, r values. 



ometric dependence. In this case the best fit corresponds 
to EPS09 parameter set 30 and <7f, r — 3.4 mb. Again, 
even this best fit has af = 46.3 which corresponds to 
an extremely poor fit. One notable feature that may 
appear counter-intuitive is that for some EPS09 nPDF 
parameter sets, the best fit shown by the red curve ap- 
pears very far below the Rcp data points. Because the 
rapidity shape is so poorly matched, it is possible that 
a better fit is obtained with the x 2 under the assump- 
tion that the global scale uncertainty of 8.2% has a three 
standard deviation fluctuation low. 



IV. NEW GEOMETRIC CONSTRAINTS 

In [17], the PHENIX collaboration presented a new 
way of using the experimental data to test and con- 
strain the geometric dependence of the combined nu- 
clear effects. By plotting the Rcp value (0-20%/60- 
88%) versus the geometry-averaged nuclear modification 
■R(jAu(0-100%), there are constrained parametric depen- 
dencies for the linear, exponential, and quadratic cases. 
In |17j . the analytic parameterization for M{r T ) as a 
function of the average A(r T ) was used to compare the 
nuclear modification in all centralities for a given param- 
eter a value. For a given geometric dependence, varying 



the values of a results in a locus of points for a con- 
strained relationship between RdAu and Rcp- Shown in 
Figure [10] as dotted, solid, and dashed black lines are 
the result of that analytic calculation for the exponen- 
tial, linear, and quadratic cases respectively. To be clear, 
these curves are calculated purely from a Monte Carlo 
Glauber geometry, the average density-weighted nuclear 
thickness A(r T ), and the simple geometric dependence 
equation (i.e. no specific model of nPDFs, ov, etc). Also 
shown are the PHENIX experimental data with the lines 
as point-to-point uncorrelated uncertainties and the el- 
lipses as one-standard deviation contours from the com- 
bined systematic uncertainties. As stated in [T7], this 
demonstrates that the forward rapidity J / tp data cannot 
be reconciled with an exponential or linear geometric de- 
pendence for the nuclear modification. 

We test this picture further by additionally plotting the 
results from our calculations using the EPS09 nPDFs and 
a br model. In Figure [lUJ we show all EPS09 nPDF pa- 
rameter sets using the exponential geometric dependence, 
a range of values of Ubr (from 0-18 mb in 2 mb steps) and 
the full range of rapidity values (as red lines) . The sub- 
set corresponding to a br = are shown as blue lines. As 
expected, since the nPDF dependence is exponential, the 
blue lines fall almost perfectly on the analytic pure ex- 
ponential case (dotted black line). With the additional 
abr contribution which also has an exponential geometric 
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FIG. 12. The points are the PHENIX J/tp R C p versus R dA u- 
The ellipses are the one-standard deviation contours for the 
systematic uncertainties. The open circles, closed squares, 
and closed circles are from backward, mid, and forward ra- 
pidity respectively. The black curves are analytic calculations 
assuming a purely exponential, linear, and quadratic geomet- 
ric dependence for the nuclear modification 17 . The blue 
lines are the full calculation, using a quadratic thickness de- 
pendence for the shadowing, for all EPS09 nPDF parameter 
sets with Obr = 0. The red lines are the full calculation for all 
EPS09 nPDF parameter sets with all 0}, r values. 



dependence, we expected that everything would collapse 
onto the same line. However, with two competing effects 
-RdAu(0-100%) =1 does not always equate with the trivial 
case of no modification, but can also occur if the two ef- 
fects average to 1. In the latter case, Rcp need not be 1. 
Specifically in our case, in the backward rapidity region 
the nPDF leads to an enhancement (anti-shadowing) and 
the Obr to a suppression. This competition can lead to 
the case where i?dAu(0-100%) = 1, while the Rcp ^ 1 
(ie. some slight enhancement in peripheral events due to 
the nPDF effect and some slight suppression in central 
events due to the Obr effect). This effect leads to the 
slight splitting of the lines for values near RdAu = 1- 



Shown in Figure 11 are the same quantities for the 
nPDF linear case in our calculation. Again, the case 
with Obr = leaves only the purely linear nPDF and thus 
the blue lines collapse onto the analytic linear case (solid 
black curve). The red curves for all cr& r > cases result 
in a geometric dependence that is part linear and part 
exponential. Thus, one sees that for larger suppressions 
(also larger Ob r values) , the red curves move between the 
analytic linear case to the analytic exponential case. One 
again sees some cases of i?dAu(0-100%) = 1, while the 
Rcp is not equal to one for the same reason as described 
above. 

Lastly, in Figure [12] we show the quadratic case. In 
this case with o\, r = the blue lines are close to the black 



dashed analytic quadratic case, but not perfectly. This 
is due to the inclusion of fluctuations in the thickness 
included in our full calculation shown here. As one in- 
creases the value of Obr in 2 mb increments one sees the 
red lines moving to the left as the exponential geometric 
dependence from Obr dominates over the quadratic nPDF 
effect. 

This full suite of curves reveals that even attempting 
to fit just the forward rapidity data with a larger and 
larger Ob r will not successfully capture the full centrality 
dependence (even using the quadratic nPDF contribu- 
tion) . 

The results of fitting Rqp versus rapidity shown in Fig- 
ures |8] and [9] demonstrated that no variation in the model 
(e.g. EPS09 nPDF parameter sets, nPDF geometric de- 
pendence, single Obr values, etc.) can be reconciled with 
the full rapidity and centrality dependence of the exper- 
imental data. It is possible that Obr may have a rapidity 
dependence due to the different relative velocity of the cc 
pair with respect to the target nucleons. A naive expec- 
tation is that a shorter time spent in the nucleus should 
result in a smaller Obr due to a smaller growth in the phys- 
ical size of the cc pair towards that of the final state J/ip. 
However, in order to attempt to reconcile the data with 
a rapidity-dependent Obr the cross section would need 
to be much larger at larger rapidity (as shown in |25] 
in Section 5.2). In fact, the results from the Rcp and 
RdAu correlation indicate that no value of Ob r (even with 
an independent value for forward rapidity) can properly 
describe the centrality dependence. This fact can be un- 
derstood by noting that no matter how large a value of 
Ob r one uses, the geometric dependence is always expo- 
nential. 

Another possibility is that the geometric dependence 
of the nPDF varies with rapidity. While this is possible, 
without direct theoretical guidance or other experimen- 
tal constraints, one would just be arbitrarily fitting the 
data. In the next section we consider two additional mod- 
els in the literature that have different physics than this 
framework. 



V. ADDITIONAL MODEL COMPARISONS 

One proposal for physics beyond the factorized nPDF 
and Obr involves calculations incorporating gluon satu- 
ration (non-linear evolution of the gluon distributions at 
low x). In the PHENIX paper [17], the data were com- 
pared with a color glass condensate calculation [SB] that 
incorporated suppression at low x from gluon saturation 
and enhancement from double-gluon exchange diagrams. 
More recent calculations following this framework [26] 
include a more accurate treatment of the nuclear geom- 
etry and the dipole-nucleus scattering amplitudes, and 
are consistent with recent calculations for the Au+Au 
1291 . Shown in Figure 13 (left panel) is the calcu- 



case 



lated J/ip nuclear modification for different rapidities as 
a function of the d+ Au event impact parameter (b) . It is 
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FIG. 13. (Left Panel) Shown are J/tjj nuclear modification results from the color glass condensate calculation including gluon 
saturation from 26 as a function of event impact parameter in d+Au reactions. (Right Panel) Shown are the results of a 
coherence and color transparency calculation [27] for the J/ip nuclear modification as a function of event impact parameter in 
p+Au reactions. 



notable that for large b there is a 30% nuclear enhance- 
ment, even though the coherence is only an effect in the 
longitudinal direction and the local nuclear density for 
these large-& events is small. 

In an alternative calculation presented in |27| , the J/tp 
production is controlled by coherence and color trans- 
parency effects. Shown in Figure 13 (right panel) is the 
calculated J/tp nuclear modification factor as a function 
of event impact parameter in p+Au reactions. 

In either case, we can fold these dependencies with 
the rf+Au event impact parameter distributions or the 
r T distribution (for the p+Au predictions), and compute 
the RdAu and i?cp modification factors to compare with 
the PHENIX experimental data. These comparisons are 
shown in Figure |14| We also show for comparison the 
EPS09 nPDF default with linear geometric dependence 
and Obr = 4 mb. The calculation from [2 7) yields simi- 
lar results to the EPS09 nPDF and o\, r calculation, and 
has insufficient suppression at the most forward rapid- 
ity. The color glass condensate calculation shows better 
agreement at forward rapidity, though the rapidity de- 
pendence is not as steep as that of the Rep experimen- 
tal data. But it has substantial disagreement with the 
data at mid and backward rapidities. Note that that cal- 
culation assumes coherence over the entire longitudinal 
extent of the nucleus, but this coherence approximation 
is no longer valid at some higher-a; values (ie. moving 
towards mid and then backward rapidity). It is also no 
longer valid for low densities that occur at large impact 
parameter values. Thus, the result that for large im- 
pact parameters the color glass condensate calculation 
approaches RdAu ~ 1.3 (shown in the left panel of Fig- 



ure |13|) is likely to simply be outside the range of validity 
of the calculation. 



VI. INITIAL-STATE ENERGY LOSS 

Another physical effect that has been proposed is that 
of initial state parton energy loss. The very forward ra- 
pidity J/tj) are produced from a high x\ parton (from the 
deuteron) and a low xi parton (from the gold nucleus). 
If the high x\ parton from the deuteron loses some en- 
ergy before the hard scattering, this will result in a lower 
J/tp production probability and a shift backward in ra- 
pidity for any produced particles (including the J/ip). 
This framework has been used in an attempt to recon- 
cile lower energy Drell-Yan data in p+A collisions (see 
for example [301 EI])- However, the same data have also 
been interpreted in terms of nuclear shadowing models 
instead of large initial-state energy loss. It is unclear if 
these experimental data sets can be reconciled with the 
same energy loss, and whether that provides a consistent 
picture for all such p(d)+A data. 

More recently, in [32 a calculation is presented of 
initial-state parton energy loss and its impact on Drell- 
Yan production with predictions for measurements in 
p + A collisions. In the case of initial-state radiative en- 
ergy loss, they predict that AE/E cx L, where L is the 
path through the nucleus prior to the hard scattering. 
Data from experiment E906 at Fermilab will directly ad- 
dress this prediction, and in the clean Drell-Yan channel 
without final state effects [3"3"j . 

As a preliminary investigation of the impact of initial- 



10 




Centrality 60-88% 
ISIobal Scale Uncertainty ±10.0% 



CGC Calc. (Tuchin) 
Coherence Calc. (KopelioviclY 
EPS09 + sigma =4 mb 



-h 




_Centrality 0-20% 

IGlobal Scale Uncertainty +8.5% 



ffl - 




Centrality 0-20%/60-88% 
'_ Global Scale Uncertainty +8.2% 



FIG. 14. Shown are the PHENIX experimental data for 
J/V> BiiAu peripheral (top) RdAu central (middle), and Rcp 
0-20%/60-88% (bottom) as a function of rapidity. The green 
curve results from [27] with coherence effects and color trans- 
parency. The blue curve is the color glass condensate calcula- 
tion [25] . The red curve is our calculation with nPDF EPS09 
default set=l with a linear geometric dependence and otr = 4 
mb. 



state parton energy loss, we have implemented this mech- 
anism in our calculation (in addition to the nPDF and 
Obr contributions). Within the Monte Carlo Glauber, 
we also calculate N tu be[before], the number of Au nu- 
cleus nucleons in the tube that have a z location prior 
to the z position of the binary collision of interest. We 
posit that the initial-state energy loss is proportional to 
Ntube [before] , which is the same as being proportional to 
the path L, with the inclusion of local fluctuations. For 
this calculation, we have utilized the PYTHIA produc- 
tion g + g — >• J/ip + X kinematics. For each binary colli- 
sion location, we randomly select a PYTHIA x±, X2 com- 
bination and the average J/tj) rapidity for those kinemat- 
ics. We then calculate the expected x\ shift due to the 
energy loss corresponding to the particular N tu be [before] 
value. From this information, we calculate the decrease 
in the probability for these partons to produce a cc pair 
and the new (lower) average final state J/tj) rapidity with 
the modified parton kinematics. We have varied the pro- 
portionality constant for the initial-state energy loss in 
the N tu be[before] dependence. 

In Figure [15] (left panel) , we show the calculation re- 
sults including only initial state parton energy loss (i.e. 
no nPDF modification and ov = 0). One observes a 



larger suppression at foward rapidity, and in fact a mod- 
est enhancement at backward rapidity. In Figure 15 
(right panel), we show the results when now assuming 
a quadratic path dependence for the energy loss (i.e. 
AE/E oc L 2 ). In this case, for large coefficients in the 
proportionality constant, there is large suppression at all 
rapidities. For either the L or L 2 dependence, one cannot 
achieve good agreement with the experimental data with 
initial-state parton energy loss alone. 

We then include all three nuclear effects (nPDF mod- 
ification, <7f, r , and initial-state parton energy loss), and 
vary the EPS09 nPDF parameter set and fit for the best 
Obr and energy loss coefficient. In Figure 16 (left panel), 
we have found the best x 2 from the fit to the Rcp- The 
blue curve represents the best overall fit for all EPS09 
parameterizations and the optimal ov and initial-state 
energy loss coefficient. The best fit gives a reasonable 
description of the Rqp and corresponds to EPS09 pa- 
rameter set 23, Obr = 3 mb, and AE/E w 0.05/fmxL 
(converting the average N tu be [before] to a length through 
normal nuclear matter density). The x 2 = 20.2 which is 
a better fit than without the initial-state energy loss, but 
still gives a probability of less than 5%. However, more 
striking is that with a reasonable fit to the Rqp, there 
is no global agreement with the overall suppression for 
RdAu in peripheral or central events. 



In Figure 16 (right panel), we show the same quanti- 



ties, but now assuming that the initial-state energy loss 
is quadratic in the path or N tu be[before]. In this case 
the best fit to Rqp has EPS09 parameter set 5, ov = 4 
mb, and the initial-state energy loss corresponds to ap- 
proximately AE/E sa 0.005/fm 2 x L 2 . Again there is no 
overall global good fit for RdAu and Rcp- 

These results are a first look in comparing the sim- 
plest initial-state energy loss calculation to these J/tp 
data. We have not included Poisson fluctuations of the 
radiated quanta, which are important when one is near 
the very high- a; i limit as pointed out in (35] ■ However, 
in their calculation they have not included the important 
fluctuations in the L value itself, as we have done through 
utlizing N tu be[before]. It is premature to draw any firm 
conclusions about the implications for initial-state par- 
ton energy loss from these comparisons with J/tp data 
alone. Further calculations are necessary, as well as com- 
parisons to observables from other collision energies and 
other final state produced particles. 



VII. SUMMARY 

In this paper, we have presented an extension of calcu- 
lations for J/tp nuclear modifications including modified 
parton distribution functions (nPDFs) and fit parameter 
a br - Utilizing the full set of EPS09 nPDFs and three dif- 
ferent postulated geometric dependencies, we find that 
the calculations cannot be reconciled with the full ra- 
pidity and centrality dependence of the PHENIX d+Au 
J/tp data for any nPDF variation and any ov value. Ad- 
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FIG. 15. (Left Panel) Calculation including initial-state parton energy loss only and with AE/E oc L or equivalently in 
our case with fluctuations to Ntube [before] . The curves correspond to coefficients of 0.01/fm, 0.03/fm, 0.05/fm, 0.07/fm, and 
0.09/fm (from upper to lower in order). (Right Panel) Calculation including initial-state parton energy loss only and with 
AE/E oc L 2 or equivalently in our case with fluctuations to N tu be[before] 2 . The curves correspond to coefficients of 0.005/fm 2 , 
0.015/fm 2 , 0.025/fm 2 , 0.035/fm 2 , and 0.045/fm 2 (from upper to lower in order). 



ditionally, the comparison with Rqp versus RdAu indi- 
cates that even a much larger aj, r at forward rapidity 
cannot reconcile the calculation with the data, since the 
o~br contribution always has an exponential geometric de- 
pendence. We also compare to two different coherence 
calculations and find no agreement across all rapidities 
and centralities. Most likely new physics beyond these 
calculations is a signficant contributor, perhaps initial- 
state parton energy loss. A first look at a simple pa- 
rameterization of initial-state energy loss allows a better 
description of the J/ip Rcp , but without a good simulta- 
neous description of RdAu- Additional constraints from 
Drell-Yan and direct photon observables at forward ra- 
pidity and at different ^/saT/v energies may be necessary 
to disentangle these effects. 
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